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I. INTRODUCTION 



The idea of combining a stochastic process for waiting times between two consecutive 
events and another stochastic process which associates a reward or a claim to each event 

'he Cramer-Lundberg 



ruling the random 



dates back at least to the first half of the XXth Century 
model for insurance risk is based on a point (or renewal) process |3, 
times at which random claims have to be paid by the company due to the occurence of 
accidents. Similar concepts have been used in renewal theory and in queueing theory, as 
well BflflSflB- 

In the 1960s, Montroll and Weiss published a celebrated series of papers on random 
walks, where they applied the ideas developed by mathematicians working in probability 
theory to the physics of diffusion processes on lattices. In particular, they wrote a paper 
on continuous-time random walks (CTKWs) |9j, in which the waiting-time between two 
consecutive jumps of a diffusing particle is a real positive stochastic variable. 

The paper of Montroll and Weiss on CTRWs was the starting point for several devel- 
opments on the physical theory of diffusion. In more recent times, CTRWs were applied 
Dack to economics and finance by Rudolf Hilfer 10], by the authors of the present paper 

■"in nn 

131 Il4j and, later, by Weiss and co-workers [15J, [16J] and by Kutner et al. 



11 



12 



17, Ui 



However, here, the focus will be on anomalous relaxation properties of the waiting-time 
probability density and on the consequent relation between CTRWs and fractional diffu- 
sion. 

Anomalous relaxation with power-law tails of the waiting-time density was investigated by 



12 1. Shlesinger, Tunaleyand other 



21 



22, 21 



sec 



means of Monte Carlo simulation by Montroll and Scher 
authors studied the asymptotic behaviour of CTRWs for large times [20 
also [25]). Hilfer has reco gniz ed the important role played by Mittag-Leffier type functions 
in anomalous relaxation ^(j . ^27 \. Interesting contributions about the theory of CTRWs 



including the problem of anomalous relaxation, can be found in refs. 



0, 0, 0, 3 



H fl 0, H 



411 ]. The recent book by ben-Avraham and Havlin discusses 



in depth the possible applications of the formalism developed in the aforementioned papers 

0. 

The asymptotic relation between properly scaled CTRWs with power-law waiting times 
and fractional diffusion processes has already been rigorously studied by Balakrishnan 
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in 1985, dealing with anomalous diffusion in one dimension |43j, four years before the 
publication of the fundamental paper b y S chneider and Wyss on the analytic theory of 



Pactional diffusion and wave equations 



45 
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44J . Later, many authors discussed this relation 



52j, |53j. As written above, the correspondence between CTRWs 



with Mittag-Leffier waiting time and the time-fractional diffusion equation has been lucidly 
worked out and explained in ref. 48j by Hilfer and Anton, who have shed light on the 
relevance of the Mittag-Leffier function, their specific aims, methods and interpretations 
being completely different from those of Balakrishnan. However, it must be recognized that 
already Balakrishnan in his formula (27) has found, as the natural choice for the waiting- 
time in CTRWs approximating fractional diffusion, the waiting-time density whose Laplace 
transform is (in the notation used in this paper) 1/(1 + cs^), where c is a positive constant. 
Implicitly, this is the Mittag-Leffier waiting-time described in Sec. Ill below. Meerschaert 
et al. have developed a method to derive the equations for CTRWs in the diffusive limit 
In their paper, they discuss both the coupled and uncoupled case. 
The present paper is devoted to a detailed discussion of the uncoupled case and it is 
organized as follows. In Sec. II, the basic quantities are introduced and a summary of 
the theory is given. Sec. Ill is devoted to the solution of the master equation in the 
uncoupled case. General formulae are presented and specialized to the case of the Mittag- 
Leffier waiting-time survival probability, in which an exact solution is available in terms of 
a fractionally generalized compound-Poisson process. In this section, a fractional relaxation 
equation satisfied by the Mittag-Leffier function is discussed. In Sec. IV, the proper scaling 
leading to the fractional diffusion equation is presented. The main result of this section is 
that the solution of the CTRW master equation weakly converges to the solution of a Cauchy 
problem for the fractional diffusion equation. This result is a version of the Central Limit 
Theorem and the steps for a rigorous proof are sketched. Finally, in Sec. V, the reader can 
find a discussion of some objections which are usually raised when dealing with fractional 
diffusion. Unnecessary mathematical difficulties have been avoided throughout the paper. 

II. BASIC DEFINITIONS 

As mentioned in the introduction, CTRWs are essentially point processes with reward. 
The point process is characterized by a sequence of independent identically distributed 
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(i.i.d.) positive random variables Tj, which can be interpreted as waiting times between two 
consecutive events: 

n 

tn = to + S ^T i '- l t n -t n -i = T n ] n = 1,2,3,...; t = 0. (1) 

i=i 

The rewards are (i.i.d.) not necessarily positive random variables: £j. In the usual physical 
intepretation, the represent the jumps of a diffusing particle (the walker), and they can 
be n-dimensional vectors. In this paper, only the 1-dimensional case is studied, but the 
extension of many results to the n-dimensional case is straightforward. The position x of 
the walker at time t is (with N(t) = max{n : t n < t} and x(0) = 0): 

N(t) 

s(f) = 5>. (2) 
i=i 

CTRWs are rather good and general phenomenological models for diffusion, including 
anomalous diffusion, provided that the time of residence of the walker is much greater 
than the time it takes to make a jump. In fact, in the formalism, jumps are instantaneous. 

In general, the jumps and the waiting times are not independent from each other. Then, 
the random walk can be described by the joint probability density </?(£, r) of jumps and 
waiting times; </?(£, r) d^ dr is the probability of a jump to be in the interval (£, £ + <i£) and 
of a waiting time to be in the interval (r, r + dr). The following integral equation gives the 
probability density, p(x,t), for the walker being in position x at time t, conditioned by the 
fact that it was in position x = at time t = 0: 

rt r+co 

p(x,t) = 5(x)V(t) + / <p(x - x',t -t')p(x',t')dt' dx', (3) 

JO J~oo 

where ^(r) is the so-called survival function, ^(r) is related to the marginal waiting-time 
probability density ip{r). The two marginal densities ipij) and A(£) are: 

/+oo 
-oo 

/•OO 

A(0 = / ^,r)dr, (4) 
Jo 

and the survival function ^(r) is: 

^(r')dr' = J ^(T')dr'. (5) 
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The integral equation, eq. can be solved in the Laplace- Fourier domain. The Laplace 
transform, g(s) of a (generalized) function g(t) is defined as: 

f +00 



= / dte-tgit), (6) 
./o 

whereas the Fourier transform of a (generalized) function /(x) is defined as: 

r+00 

/(*) = / dxe^f(x). (7) 



A generalized function is a distribution (like Dirac's 5) in the sense of S. L. Sobolev and L. 

n 

Schwartz j54| . 
One gets: 



p(K,s) = *(s) i- -, (8) 

1 - ip{K,S) 



or, in terms of the density ip(r): 



as, from eq. (jSJ, one has: 



'-* W ' . (9) 

s l-(p{K,S) 



9{s) = (10) 
s 



In order to obtain p(x, t), it is then necessary to invert its Laplace- Fourier transform p(n, s). 
Analytic solutions are quite important, as they provide a benchmark for testing numer- 
ical inversion methods. In the next section, an explicit analytic solution for a class of 
continuous-time random walks with anomalous relaxation behaviour will be presented. It 
will be necessary to restrict oneself to the uncoupled case, in which jumps and waiting-times 
are not correlated. 



III. SOLUTION OF THE MASTER EQUATION 

In this section, the solution of eq. (J3J) will be discussed in the uncoupled case. First of 
all, a general formula will be derived for p(x, t), then it will be specialized to two cases: the 
well-known case of an exponential survival function and the case where the survival function 
is a Mittag-Lefner function. The connections and the analogies between these two cases 
will be presented. A new solution will be obtained in terms of a fractionally-generalized 
compound Poisson process. 
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As anticipated above, the study will be restricted to uncoupled continuous-time random 
walks. This means that jump sizes do not depend on waiting times and the joint probability 
density for jumps and waiting times can be factorized in terms of the two marginal densities: 



with the normalization conditions J c?£A(£) = 1 and J drijj(r) = 1. 
In this case the integral master equation for p(x, t) becomes: 

ft r ;>+oo 

p(x,t) = 5{x)^{t) + ip(t-t') / \{x-x')p(x',t')dx' 

JO \_J-oo 



(11) 



dt' 



(12) 



This equation has a well known general explicit solution in terms of P(n, t), the probability 
of n jumps occurring up to time t, and of the n-fold convolution of the jump density, A n (x): 



oo r+co 



X n (x) = / . . . / d£„_id£n-a . . . d^X(x - £ n -i)A(£„-i - £»- 2 ) • • • A(£i). (13) 

J —00 J —00 J— 00 

Indeed, P(n, t) is given by: 

P(n,t)= [ tp n (t-r)^(r)dr (14) 
J 

where ip n { T ) is the n-fold convolution of the waiting-time density: 



rfr„_idr n _ 2 • • • dnip(t - T n _i)V>(r n -i - r n _ 2 ) • • • if>{ri). (15) 



'0 Jo Jo 

The n-fold convolutions defined above are probability density functions for the sum of n 
variables. 

The Laplace transform of P(n,t), P(n,s), reads: 



P(n,s) = [^(s)] n ^(s) 

By taking the Fourier-Laplace transform of eq. (|12p. one gets: 

1 



(16) 



p(n,s) = *(s) — 

l-V(s)A(«) 

But, recalling that |A(«)| < 1 and \i/j(s)\ < 1, if k ^ and s 7^ 0, eq. (|T7|) becomes: 

00 

n=0 



(17) 



(18) 



this gives, inverting the Fourier and the Laplace transforms and taking into account eqs. 
(H3J) and (JUJ): 

oo 

p(x,t)=^ / P(n,t)X n (x) (19) 

n=0 

Eq. (JT5j) can also be used as the starting point to derive eq. (|T2*j) via the transforms of 
Fourier and Laplace, as it describes a jump process subordinated to a renewal process. 

A remarkable analytic solution is available when the waiting-time probability density 
function has the following exponential form: 

^(r) = /ie- MT . (20) 

Then, the survival probability is ^(t) = e _AtT and the probability of n jumps occurring up 
to time t is given by the Poisson distribution: 



P(n i t) = ^-e-' a . (21) 



nl 



In this case, equation (fTT?|) becomes: 

P(x,t) = jr i -^fe-» t \ n (x). (22) 

n=0 

When A (re) is the jump density for a positive random variable, eq. (J22J) is the starting point 
of the Cramer-Lundberg model for insurance risk . It is worth noting that the survival 
probability ^(t) satisfies the following relaxation ordinary differential equation: 

A^(r) = -(iV(r), r > 0; ^(0 + ) = 1. (23) 
dr 

The simplest fractional generalization of eq. ()23j) giving rise to anomalous relaxation 
and power-law tails in the waiting-time probability density can be written as follows, by 
appropriately choosing the time scale: 

dP 

_v[/( r ) = _H/( r ), r > 0, < (3 < 1; *(0+) = 1, (24) 

where the operator d^/dt 13 is the Caputo fractional derivative, related to the Riemann- 
Liouville fractional derivative. For a sufficiently well-behaved function f(t), the Caputo 
derivative is defined by the following equation, for < (3 < 1: 

d^ m = Txr^)dtj (^7F dr -f(r^) /(0 } ' (25) 



and reduces to the ordinary first derivative for (3 = 1. The Laplace transform of the Caputo 
derivative of a function fit) is: 

C (^ mS ) = s "f( s )- s "~ 1 f(U + )- ( 26 ) 
If eq. (f2T)|) is applied to the Cauchy problem of eq. (J2I|), one gets: 

*W = i^jr (") 

Eq. (|27jl can be inverted, givi ng the solution of eq. (|2*%|) in terms of the Mittag-Lemer 



function of parameter f3 
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*(r) = E p {-rfi), (28) 
defined by the following power series in the complex plane: 

oo n 

For small r, the Mittag-Leffler survival function has the same behaviour as a stretched 
exponential: 

*(r) = ^(-r") ~ 1 - ^ ~ exp{-r^/r(/3 + 1)}, < r « 1, (30) 
whereas for large r, it has the asymptotic representation: 



$ (t) ^M )0 < /3 <1 ) ^oo. (31) 

7T 7~" 



Accordingly, for small r, the probability density function of waiting times ip(r) = — d^(r) / cir 
behaves as: 

1>(T) = -frfyi-T*) ~ T (J3 + 1) expj-r^/r^ + 1)}, < r « 1, (32) 
and the asymptotic representation is: 



sm(f3ir) T(p + 1) 

^(r) ^ % +1 j , < p < 1, r - oo. (33) 



Before going on, it is now time to review the results obtained so far. The solution of 
equation (|24|) is a survival probability function \I/(t) with power- law decay r _/3 if < (3 < 1 
and r — > oo. The decay exponent of the corresponding probability density function ip(r) is 

8 



— (/3 + 1), with values in the interval (1,2). This ensures that the normalization condition 
for ip(r) can be satisfied. However, already the first moment of ip(r) is infinite. It is worth 
stressing that the case (3 = 1 does not correspond to a r -1 -decay of the survival probability, 
but to the exponential relaxation described by eq. (1231) . 
The Laplace transform of ip(r) is given by (see eq. (|1U|1): 

$(s) = 1 - a^(a) = — ^. (34) 

Therefore, recalling eq. (fT7)|) and eq. (|27J) . one can obtain the Laplace transform of P(n, t): 

1 s?- 1 

P(n,s) = — -5T— -j. (35) 



This can be analytically inverted as (see eq. (1.80) in ref. [57]): 

lis?- 1 

where: 



Clt*E£H-*);-)= (36) 



Eq. yields an explicit analytic expression for P(n,t): 

t f3n 

n\ 

Eq. (|3"7j) generalizes the Poisson distribution (|21|) for the anomalous relaxation case under 
study (0 < (3 < 1). It reduces to the Poisson distribution in the case (3 = 1, in which 
the Mittag-Leffler function coincides with the exponential function. As an immediate con- 
sequence of this result and of eq. (fTHjl. one also gets the analytic solution of the master 
equation (|12|) for a continuous-time random walk characterized by the survival function of 
eq. (|2H1): 

p(^*) = E 1 r4 n) (-^) A »( a? )- ( 38 ) 

ra=0 



P(n,t) = ^W(-0 (37) 



As a consistency check, one can show that 

/+oo 
p{x,t)dx = 1, Vt (39) 
-oo 

This is equivalent to the requirement that the Fourier transform computed in k = satisfies 
p(0, t) = 1, Vt. But p(0,t) is given by: 

t /3n 



n=0 



^^E^^-^ (40) 



and recalling that for any sufficiently well-behaved function, /: 

n=0 

identifying a = —t 13 and 5 = +i , one has the following chain of equalities: 

00 f f3n 

ft '*) = = MH p ) + tP ) = Mo) = 1. (4i) 

n=0 

It is now interesting to investigate the behaviour of the exact solution given by eq. (J3%|) in 
the so-called diffusive or hydrodynamic limit. This limit is obtained by making smaller all 
waiting times by a positive factor r, and all jumps by a positive factor h and then letting r 
and h vanish in an appropriate way. This will be the subject of the next section. 



IV. THE DIFFUSIVE LIMIT 



In this section, for the first time, a collection of results by the authors of this paper 
is made available in a complete way; mathematical subtleties have been recalled wherever 
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53 



581 . Here, the focus is on the well- 



necessary. Partial results were discussed in refs. 
scaled transition to the diffusive limit based on sound limit theorems of probability theory. 
The following derivation should help the reader in judging whether, in the problem he/she 
is dealing with, the connection between CTRWs and fractional diffusion is relevant. 

As mentioned above, in order to discuss the diffusive limit, the waiting times are multi- 
plied by a positive factor r so that one gets: 

t n (r) = m + rr 2 + . . . + rr n . (42) 

Analogously, the jumps are multiplied by a positive factor h. Letting Xo(h) = 0, one has: 

x n {h) = h£ 1 + h& + ... + h£ n . (43) 

The probability density function ?/v(r) of the scaled waiting times is related to ip(r) in the 
following way: 

M r) = r > 0, (44) 

The scaled-jump probability density function \h{£) is given by: 

MO ^ (45) 

10 



The Laplace transform of ^v(t) and the Fourier transform of A^(£) are, respectively: 

Ms) = $(rs), A h («) = A(/wc). (46) 
In the Fourier-Laplace domain, the rescaled solution of the master equation reads: 

^(/C 5) = p^- ~ ~ • (47) 

S 1 - 1p r {s)X h {K) 

Eq. (|47j) will be the starting point for investigating the diffusive limit of the solution pre- 
sented in eq. (J2EJ). The results discussed above, from eq. ()42jl to eq. 1)47)1. are rather general. 
It is now possible to specialize them to the class of waiting-time densities discussed in Sec. 
Ill and to a large class of jump densities. 

For < /3 < 1, eq. ()H3)) gives the asymptotic representation of the waiting-time density. 
For such a behaviour, one has, for each fixed s > that: 

Vv(s) = V^(rs) = 1 - cifrs)' 3 + o(r p ), r -> 0. (48) 

In the case under study, it turns out that c\ — 1. Remarkably, this result holds also for 
(3 = 1. An important class of symmetric jump densities (A(— £) = A(£)) is characterized by 
the following behaviour, for 6 > and some parameter a G (0, 2): 

A(x) = [6 + e(|x|)]|x|- (a+1) , (49) 

with e(|a;|) — > as |a;| — > oo. For these densities, exhibiting a power-law decay at infinity, 
the following asymptotic relation holds: 

A h («) = A(/ik) = 1 - c 2 (%|) Q + o(/i a ), h->0, (50) 

where the constant c 2 is given by: 

c bn (51) 

2 r(a + 1) sin(a7r/2) 1 ; 

Eq. (j5D)) is valid also for symmetric densities with finite second moment a. In that case, one 
has a = 2 and c 2 = a 2 /2. Both the results in eq. ([48)1 and in eq. (|5U|) are less trivial than 
they seem. Indeed, in order to prove eq. ()48j) . it is necessary to use a corollary on Laplace 



transforms discussed in the classical book by Widder (see ref. 
proof of eq. (1501) is tedious but can be distilled from chapter 8 of 

n 

proof can be found in ref. [58J. 
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p. 182), whereas the 
60(. A simpler but weaker 



By using the asymptotics in eq. (|4*%)) and in eq. (j3Uj) and replacing in eq. (|37jl . it follows 
that: 

rC~h(K,s)~ -t—z — — — , r, h— >0. (52) 

Now, the following scaling relation can be imposed: 

c ir P = c 2 h a , (53) 

yielding, for r, h 0: 

p r , ft (/c,a)-> 8/?+[/c|tt . (54) 

This limit coincides with the Laplace-Fourier transform of the Green function (or funda- 
mental solution) for the following fractional diffusion Cauchy problem: 

0^u(x,t) = —-—u(x,t), 0<a<2, < (3 < 1, 
u(x, + ) = 6(x), xG(-oo,+oo), t > 0, (55) 

where d^/dt^ is the Caputo derivative defined in eq. (|25jl. and 9 a /9|x| Q is the Riesz deriva- 
tive: a pseudo-differential operator with symbol — |k| q . Recalling eq. (|26)1. the Laplace- 
Fourier transform of u(x,t) reads: 

s) = — ; — — , (56) 

and therefore, as anticipated, one has, for r,h — > under the scaling relation eq. (|53j) : 

p^ h {K,s) -»• u(«,s). (57) 

In this passage to the limit, p^^K, s) and «(/«, s) are asymptotically equivalent in the Laplace- 
Fourier domain. Then, the asymptotic equivalence in the space-time domain between the 
master equation eq. (|T2*|) and the fractional diffusion equation, eq. is ensured by 

the continuity theorem for sequences of characteristic functions, after the application of the 
analogous theorem for sequences of Laplace transforms [7]. Therefore, there is convergence in 
law or weak convergence for the corresponding probability distributions and densities. Here, 
weak convergence means that the Laplace transform and/or Fourier transform (characteristic 

unction) of the probability density function are pointwise convergent (see for details ref. 

3]). In other words, under the appropriate scaling, defined by eq. and in the limit 
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r,h — > 0, the solution given in eq. 
fractional diffusion equation eq. ([51))): 



weakly converges to the Green function of the 



u(x, t) 



where W a ^{u) is given by: 



+00 



-W, 



a, 



X 



(58) 



W a , p {u) = — 



dn e 



-IKU 771 



35 



53, 



61 



(59) 

62j | . In the case j3 — 1 



the inverse Fourier transform of a Mittag-Leffler function 
and a = 2, the fractional diffusion equation reduces to the ordinary diffusion equation and 
the function H / 2,i( , u) becomes the Gaussian probability density function evolving in time 
with a variance a 2 = It. In the general case (0 < (3 < 1 and < a < 2), the function 
W a> p(u) is still a probability density evolving in time, and it belongs to the class of Fox 
H-type functions that can be expressed in terms of a Mellin-Barnes integral as shown in 
details in ref. [si | . 

The scaling equation, eq. (|53jh can be written in the following form, where C is a constant: 



h = Cr Pla . (60) 

If /3 — 1 and a = 2, one recognizes the scaling typical of Brownian motion (or the Wiener 
process). Indeed, this is the limiting stochastic process for the uncoupled continuous-time 
random walks with exponential waiting-time density and the class of jump densities with 
finite second moment. In all the other cases considered in this paper: (3 G (0, 1) and 
a G (0,2), the limiting process has a probability density function given by u(x, t) in eq. 
(HE!). 



V. DISCUSSION AND CONCLUSIONS 



In this paper, the connection between a class of CTRWs with Mittag-Leffler survival 
probability and the fractional diffusion equation has been discussed. In Sec. Ill, an explicit 
solution of the master equation has been derived for long-tail processes with Mittag-Leffler 
survival function. As shown in Sec. IV, it turns out that, for this class, the solution of 
the CTRW master equation weakly converges to the solution of a Cauchy problem for the 
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fractional diffusion equation. Although such weak convergence also occurs for the waiting- 
;ime densities whose Laplace transforms have an s — > asymptotics 1 — c\s^ + o(s^) (see refs. 



43, 



581]). the Mittag-Lefner waiting-time law deserves special attention as, without passage 



to the diffusion limit, it leads to a time- fractional master equation, just by insertion into 
the CTRW integral equation. This fact was discovered and made explicit for the first time 
in 1995 by Hilfer and Anton Therefore, this special type of waiting-time law (with its 
particular properties of being singular at zero, completely monotonic and long-tailed) may 
be best suited for approximate CTRW simulation of fractional diffusion. It must be stressed 
that both the results of Sec. Ill and IV are based on sound and original mathematical 
considerations. 

It is important to remark that the presence of the time Caputo fractional derivative (or 
equivalently of the Riemann-Liouville derivative) and of the space Riesz derivative in eq. (j53|) 
is a natural consequence of the well scaled diffusion limit discussed in Sec. IV. This should 
be already clear from previous work on the relation between CTRWs and fractional diffusion 
(see, in particular ref. [48]). However, it is still often argued that there is an arbitrariness 
in the choice of the fractional operator in eq. 1)55)1 . If one uses different fractional operators, 
the physical meaning, if any, of these versions of eq. ()55j) will be different. 

Another point has been raised on the physical meaning of eq. ((55)1 . In particular, some 
authors consider the space fractional derivative unphysical due to its non-locality. An answer 
to this objection is that it is always possible to use an equation as a phenomenological model 
if it gives good results in the description of a physical phenomenon. For instance, the usual 
Fourier diffusion equation is not invariant for time inversion, whereas the basic equations 
of classical mechanics are. Still, the Fourier equation gives very useful results when used in 
many applications. 

Finally, it is important to discuss some recent results by R. Hilfer He has 

shown that not every continuous-time random walk with a long time tail is asymptotically 
equivalent to a diffusion equation with a fractional time derivative. In 0], he considers 
different ways to define fractional derivatives in time. He shows that only the Caputo type 



leads to mass-conserving fractional diffusion. In [64J, an example of a CTRW has been 
given whose waiting-time density has a power-law behaviour but whose diffusive limit is not 
the time-fractional diffusion equation. This latter counterexample seems to be in contrast 
with what has been said in Sec. IV above. However, the paradox disappears if one takes 
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into account the proper scaling given by eqs. (j44H46j) . Indeed, the counterexample by Hilfer 
does not satisfy this scaling. More precisely, the non-relevance of this counterexample for the 
theory developed in Sec. IV can be stated as follows. The waiting-time density of the second 
model presented by Hilfer cannot be written in the form of eq. (J44j) : ip r { T ) — ip(r/r)/r. 
Essentially, each of the two addends of Hilfer's density has a different scaling form. The 
scaling of eq. (jjjj) has already been used by the present authors in ref. j3| . It was previously 
used by Feller in deriving the diffusion equation from the simple symmetric random walk 



651 ] . by Balakrishnan in 1985 |43( and, in recent times, independently from the authors of 



this paper, by Uchaikin and Saenko 

We are currently working to extend our approach to the coupled case, but this will be 
the subject of future papers. 
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